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NEWTON-RAPHSON— CONJUGATE -GRADIENT TECHNIQUE FOR 
CALCULATION OF AXISYMMETRIC BUCKLING OF SHALLOW 
SPHERICAL SHELLS WITH VARIABLE EDGE CONSTRAINT* 

By Harry G. Schaeffer 
Langley Research Center 

SUMMARY 

The feasibility of utilizing mathematical programing techniques for determining 
stable equilibrium states of nonlinear structural systems is explored. The physical prob- 
lem chosen for evaluation of this methodology was the nonlinear behavior of a shallow 
spherical shell subjected to uniform pressure. The procedure found most useful was a 
hybrid method consisting principally of a generalized Newton-Raphson technique with 
intermittent application of the conjugate -gradient search technique. The variation of 
buckling pressure with changes in edge restraint and shell geometry is found for a uni- 
formly loaded shell. 


INTRODUCTION 

In the areas of nonlinear mechanics and mathematical programing, significant 
progress has been made in the development of algorithms for minimizing a function of 
several variables. Typical of problems in nonlinear mechanics for which minimization 
procedures are particularly adaptable is the determination of stable axisymmetric equilib- 
rium states of shallow spherical shells loaded by external pressure. This physical sys- 
tem exhibits the phenomenon of snap-buckling when a critical load is reached. The prob- 
lem of a spherical shell with both clamped and simply supported edges has received a 
great deal of attention in recent years (refs. 1 to 11). References 3, 7, and 8 give 
clamped edge results, and reference 5 gives simply supported edge results. However, 
since no research has been directed toward shells with edge constraints between these 
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two extremes, it would be of interest to investigate the buckling behavior for edge condi- 
tions which span these limits. 

Variational procedures using a set of assumed functions, such as the Ritz technique 
(ref. 12), have been popular methods for obtaining approximate equilibrium states of both 
linear and nonlinear systems. In addition, several investigators (refs. 13 to 16) have 
shown that an appropriate set of finite -difference equations for linear conservative sys- 
tems can be formulated by direct minimization of the energy function with respect to each 
unknown nodal displacement. 

Nonlinear finite-difference equations are derived herein in a similar manner. The 
purpose of this report is to investigate the applicability of mathematical programing 
methods as a means of determining stable equilibrium states of nonlinear systems. The 
buckling result is an example of a difficult point to obtain. Several minimization algo- 
rithms are considered, and a hybrid method centered on the generalized Newton-Raphson 
procedure is used for calculations of buckling results for the spherical shell. New 
results are presented which show the variation of the buckling pressure of spherical 
shells with changes in edge constraint. In this study only axisymmetric buckling states 
are considered. 
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SYMBOLS 


radius of curvature of spherical segment 


radius from shell axis to edge of shell (see fig. 1) 


bending stiffness, 


Et*- 


12(1 - i/2) 


nondimensional bending stiffness parameter, — 

Z 2 

Young's modulus of elasticity 


e r ,eg,e r 0 middle surface strains 

g^ gradient of the potential function II at x^ 

H maximum shell rise 


assumed location of minimum of II(x) 
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k 

k r ,kg 


l 

l 

N 

n 

P 

P 


Pci 

Per 


qi^2^3 

r 

s 

t 

u 

u 

v i 

w 

w 


nondimensional membrane stiffness parameter, 12 
middle surface bending distortions 
reference length 

nondimensional reference length, 1 / a 

number of real shell stations along a meridian (fig. 2) 

number of unknowns 

dimensional surface pressure 

nondimensional surface pressure, p/p c i 

classical buckling pressure defined in equation (5) 

axisymmetric buckling pressure corresponding to "top-of-the-knee" 
(point U in fig. 3) 

boundary restraint coefficients (spring rates) 
radius to shell middle surface from shell axis (see fig. 1) 
scalar, optimum distance to proceed in direction Vj 
shell thickness 

displacement along shell meridian 
nondimensional meridional displacement, u/Zb 
best direction to proceed from x® (eq. (11)) 
displacement normal to shell middle surface 
nondimensional normal displacement, w/Zb 
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X 


set of all Xj 


ij 


set of all nodal displacements Uj,Wj 
interval between shell stations, b/N 


preassigned small number (inequality (18)) 


Kronecker delta, 6ij = 0 for i + j; 6^j = 1 for i = j 


e i’ e i±(l/2) i n t e g ra tion factors evaluated at station i, and midway between station i 
and i ± 1, respectively 


€ •« ,€ 


r> c 0 


M 


strains 

circumferential coordinate 

geometric parameter for spherical shell, 

Poisson’s ratio 


jl2(l - t/2)b 2 
at 


n,n e ,n s ,n u 


n 


nondimensional total potential energy, edge force energy, surface 
force energy, strain energy, respectively 

dimensional potential energy 


ni,n n i nondimensional linear and nonlinear components of II, respectively 


P 


nondimensional radial coordinate, r/b 
nondimensional surface pressure, 


a>i,a >2 nondimensional curvatures; for spherical caps = Wg = ^ 
Superscripts: 

ct,{S,y denote sequence number in iterative procedures 

Primes indicate differentiation with respect to the nondimensional radial 
coordinate p. 
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POTENTIAL ENERGY FOR A SHALLOW SPHERICAL SHELL 


The shell geometry is shown in figure 1. The shape is taken to be a shallow 
spherical cap. The shell has a thickness t and a spherical radius of curvature a. 

The displacements u and w are shown positive. The displacements u and w and 
the slope 9w/ 9r may be elastically restrained at the shell boundary. The loading is an 
axisymmetric surface pressure p normal to the shell midsurface and positive outward. 
This pressure may vary along the meridian. 

The total potential energy of the shell II is composed of the strain energy n u , 
the energy of the external surface pressure n s , and the energy of boundary con- 
straints n e : 

n = n u + n s + n e (i) 

In equation (1) the potential energy is nondimensional and a bar indicates a dimensional 
quantity where 


The strain energy for a shallow spherical cap is (ref. 17) 


"u ■ £[“(' 


’22 

+ 2z^6j,00 


i) + d ( k r + k e + 2 ^ k r k 0 ) P dp 


where the strains and bending distortions are 


1 o 

e r = WjW + u' + -(w*)* 5 


e 8 = Vp 

k r = -w" 


and where the nondimensional radius, displacements, curvatures, and stiffnesses are 
given by 


5 



b 

W 1 = “2 = a 



The potential energy for a uniform surface pressure p on the shell is 

n s = -2 f fiwp dp (4) 

J 0 


where 



and 


P = 


Pci 


The classical buckling pressure p cl for a complete spherical shell is 

4Et 


p cl = 


a< 


2 \|l2 (l - v 2 ) 


(5) 


The potential energy for the arbitrary edge constraints (linear springs) integrated 
about the circumference at p = 1 is 


n e = qi u b + ^2 w b + q 3( w b)" 


( 6 ) 


where q 1? q 2 , and q 3 are meridional, normal, and rotational nondimensional spring 
constants at the boundary. The subscript b denotes the subscripted quantity at the 
boundary p = 1. 

The total potential energy from equations (1), (2), (4), and (6) is 


n = u^w + u* + |(W) 2 ] 2 + ( W 2 W " ^ + 2y [ w i w + u ' + ^ w ’) 2 ] ( W 2 W ‘ 


+ d 


(w”) 2 + (^) 2 + 2vW " W - - 20w dp + q,u 2 + q 2 w 2 + q^w^ 


( 7 ) 


It is necessary to approximate equation (7) by a finite number of unknowns. This 
is accomplished in the present investigation by approximating the continuous functions 
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(u,w) and the shell properties at N discrete stations along the meridian of the shell. 
These discrete stations are shown in figure 2. The radial interval between stations A 
is taken to be a constant, and the off-shell station N + 1 is added to allow the evalua- 
tion of derivatives at station N. 

For convenience the numerical approximation of the potential function for the 
shallow spherical cap is written as 


n = IIj + IIjji 


( 8 ) 


where 11^ involves terms in (u,w) of second degree and lower and n n i involves terms 
greater than second degree. Approximation of definite integrals by trapezoidal summa- 
tions and the use of consistent numerical approximations for derivatives give the fol- 
lowing expressions for and n n ^: 


n l = 



, 2l/U i-(l/2) U i-(l/2) 
P i-(l/2) 


+ d i-(l/2) n 


[ W, i-(1/2)| : 


p f Pi-(l/2) e i-(l/2) 

P i-(l/2) 


+ <k* 


2/22 \ 
w i V W 1 + w 2 + 2vCj0 i (Jt) 2j 


-f 

+ —n + 2^ C0«] 

Pi 


W.U. 


— 

AvjwVX 

1 

\ 

(“2 + w l) p. 

+di 

(w I') 2 + 

( v)J 

- 2S2iW. 

\ 

H 


+ 


«2 w N +q 3(*N) 


2 


(9) 


”..=1 

An advantage of the potential-energy approach is that only geometric conditions are 
required at the pole, whereas procedures based on approximating the differential equa- 
tions by finite differences require a shear condition also. These geometric conditions 
are imposed prior to minimizing the potential-energy function and require only that clo- 
sure conditions specified at the pole insure finiteness of strain. Thus, the expressions 


M»i /j / 2 \ i 

( W 1 + w 2) w i-(l/2) + u i-(l/2) + Pi . (1/2 ) > Pi-(l/2) e i-(l/2) 


( 10 ) 
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for the circumferential midplane strain eg and the midplane bending curvature kg 
(eqs. (3)) yield the conditions at p = 0, u(0) = 0, and w'(0) = 0. 

FUNCTION MINIMIZATION 
Computational Procedure 

At this point it is necessary to set up an algorithm to perform the minimization of 
the energy (as given by eq. (8)), which may be identified as the object function n. A 
computational procedure was chosen here on the basis that it might be applied to some 
general functional for many classes of problems involving a large number of variables 
where round-off error might be an important consideration and where maximums and 
inflection points may exist near the desired minimum. 

Sequential search methods of minimization (refs. 18 to 27) were considered because 
such methods are quadratically convergent. The Newton-Raphson procedure, which 
utilizes both first and second derivatives, was considered and would have been used 
exclusively except that this method cannot be guaranteed always to tend toward a minimum. 
This deficiency will be discussed further in a subsequent section. The steepest-descent 
(ref. 20), variable-metric (refs. 24 and 25), and conjugate -gradient methods all appeared 
advantageous since only first derivatives of the object function were required. The 
steepest-descent method was ruled out since it is known to exhibit poor convergence char- 
acteristics for certain classes of problems (ref. 20). The variable-metric method 
requires storage and manipulation of a matrix and therefore is not suited to problems 
involving a large number of variables. The conjugate-gradient technique was the most 
advantageous method (after the Newton-Raphson procedure) but could not be used exclu- 
sively since the convergence characteristics of this method appear to be adversely affected 
by round-off error (see, also, refs. 21 and 26). 

The linear axisymmetric behavior of a uniformly loaded circular plate was obtained 
in a sample calculation in which the plate was approximated by 22 unknowns. Since the 
potential function for this problem is quadratic, the conjugate -gradient method should con- 
verge in 22 iterations. For the plate example, however, approximately 400 iterations were 
required for single -precision calculations (8 digits) and 40 iterations for double-precision 
calculations (16 digits). This example shows the strong influence of round-off error on the 
results and ruled out the pure conjugate-gradient method for use in minimization of the 
potential function. 

Thus, the minimization technique finally chosen was that of utilizing primarily the 
Newton-Raphson method but incorporating the conjugate -gradient method when the Newton- 
Raphson method failed to move toward a minimum point. This method, which will now be 
described in some detail, was used in this investigation to determine the buckling pressure 
of shallow spherical shells subjected to external pressure. 
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In both the Newton-Raphson and the conjugate -gradient methods, a new approxima- 
tion of a dependent variable x is given in terms of the present approximation in the 
following form: 


( 2+1 a a a 

X i = X^ + S Vj 


( 11 ) 


where the superscript a denotes the iteration number. 

The vector v? represents the best direction to proceed from x“ , where x“ 
represents a point in an n-dimensional space and n is the number of unknowns. The 
scalar s a represents an optimum distance to be traversed from x“ in the direc- 
tion vf\ Here x, is the set of all nodal variables associated with u and w. Equa- 

11 1 o rv fv-j-l 

tion (11) is to be used to generate a sequence of points x£ , Xj , . . ., x“ , Xj sub- 
ject to the condition that 

n(x Q!+1 ) < n( x Q! ) (12) 


Newton-Raphson Procedure 

The algorithm selected to perform the minimization was a generalized Newton- 
Raphson procedure (ref. 18). This method utilizes both first and second derivatives of 
the object function II and assumes that II(x) and the first two derivatives are con- 
tinuous. The use of this method for the expansion of II(x) at the point xj = x“ leads 
to the following expression where terms higher than the second order are neglected: 


n(x) = n(x«) + 8n ^— ( x i 





1 a 2 n(x Q! ) 

2 9xj Qxj 




(13) 


Since a necessary condition for a relative minimum is that the gradient of n(x) vanish, 
it follows that 

an(x a ) i^nGc 01 ) a _ 0 

9xi 9xj 9xj j 


where terms greater than the first degree in Xi have been dropped and where 

a ,a a 

e i =h i ' x i 




/V QTT II /y 

x. . Since the quanties -2H- and — — — — are known at x. , the error vector e“ can 

] axi ax, ax, ] ’ ] 


axi axj 

be foimd by the matrix inverse as follows: 



(14) 


(15) 

h j 

and 

a 

e j 

can 


a 

e i=- 


a 2 n(x“) 


axi axj 


-l 


8n(x a ) 

ax« 


(16) 
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The derivative terms on the right-hand side of equation (16) are presented for the object 
function II as given in the appendix. For this problem the matrix to be inverted on the 
right-hand side of equation (16) is symmetric and has a narrow band of nonzero terms. 
The error vector e® is therefore obtained conveniently from equation (14) using 
Gaussian elimination. 


Once the error vector 


e® is calculated, then the new approximation becomes 

a+l a a 

x:~ = xf + er^ 

J J J 


(17) 


This sequence of calculations (eqs. (16) and (17)) is repeated until 


a 

ICO 

V 

X® 

] 

1 



(18) 


where 6 is a preassigned small number. When inequality (18) is satisfied, then the 
procedure has presumably converged to the solution hj where 


The generalized Newton-Raphson procedure has a drawback in that inequality (12) 
is not necessarily satisfied. The solution may tend toward a maximum or stationary 
point rather than a minimum. A check is made after each iteration to insure that inequal- 
ity (12) has been satisfied. When the test fails and inequality (12) is not satisfied, then a 
new starting point which guarantees a move toward a minimum is required. A conjugate- 
gradient technique was selected for the restart toward a minimum. 


Restart by Conjugate-Gradient Method 


Assume that at some iteration number (0 + 1) the inequality (12) is not satisfied and 
that is the last approximation satisfying the inequality (12). By the conjugate- 
gradient method (refs. 20 to 23) the unknown vector v£, required in equation (11), is 
given by 



(19) 


The distance 


s / 3 


in the direction 


is given by 


^(xf + s ^vf) = 0 (20) 

ds p ' ' 

From equation (11), 

x^ +1 = x^ + s^v^ (21) 
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For subsequent iterations by this method, 


where 


v[ = 


-g( + 


g 


g( 




g 


t -r 1 

y-1 i 




( 22 ) 


» and s^ is given by equation (20) with /3 replaced by y. Each iteration leads to a new 
approximation from equation (11). 

Only a few cycles are usually required to satisfy inequality (12), and then the 
Newton-Raphson procedure can be resumed. 


APPLICATION TO BUCKLING OF SPHERICAL SHELL 


The pressure-deflection curve for a uniformly loaded shallow shell has the general 
form shown in figure 3 (see, e.g., ref. 8). The buckling pressure p cr is defined as the 
pressure corresponding to point U on the curve. This curve OULN is representative 
of those shells for which > ll, where y? is defined as 

M 2 = \/l2(l - v 2 )(^j 

Very shallow shells (/i 2 < ll) are of no interest since the curve OULN increases mono- 
tonically and there is no definable top-of-the-knee. 

In the prebuckled range, the Newton-Raphson procedure converges rapidly when 
the linear solution is used as the initial-state point. In the region of interest for 
buckling pressures, that is, near point U, convergence difficulties may occur when the 
Newton-Raphson procedure alone is used and when the most recent converged- state point 
or the linear solution is used as an initial state. For pressures greater than that asso- 
ciated with point U with a similar assumed initial state, the Newton-Raphson procedure 
will not converge. In the region of interest approaching point U a restart by the conjugate- 
gradient method needs to be applied. 

Two convergence tests are made in the application of the Newton-Raphson procedure. 
The first (inequality (12)) determines that the function is being minimized and the second 
(inequality (18)) determines the location of the minimum to an acceptable degree of accu- 
racy. In this study 6 was taken as 0.05. 
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If inequality (12) is not satisfied after a prescribed number of restarts of the 
Newton-Raphson — conjugate -gradient procedure, the pressure increment is reduced and 
the procedure is again applied. If two successive reductions of the pressure increment 
do not lead to convergence, by satisfying inequality (12), then it is assumed that the 
pressure associated with the last converged result is the buckling pressure. In the pres- 
ent investigation the new pressure increment is taken as 1/5 of the current pressure 
increment so that the buckling pressure is found to within 0.025Ap, where Ap is the 
magnitude of the pressure increment. When the Newton-Raphson procedure is reini- 
tiated after reducing the pressure increment, the last converged solution is used as an 
initial approximation of the location of the minimum. 

EFFECT OF EDGE CONSTRAINTS ON THE AXISYMMETRIC 
BUCKLING OF A SPHERICAL SHELL 

The minimization procedure and resulting computer program were used to study 
the influence of arbitrary edge constraints on the axisymmetric buckling pressure of a 
spherical cap. Comparisons are made between present results and those of previous 
investigators. For the present investigation, 10 finite-difference increments were used 
along a shell meridian together with the following parameter values: 

v* = 0.1 
b/t = 100 
q 2 = 10* 

where the value of qg is essentially large enough to constrain normal edge displacement. 

The variation of the axisymmetric buckling pressure with shell parameter ju for 
clamped and simply supported edges was obtained for comparison with results presented 
by other investigators as an evaluation of the minimization procedure. (For example, 
see ref. 8 for clamped results and ref. 28 for simply supported results.) The results 
are shown in figure 4 where it is seen that buckling pressures as determined by the pro- 
cedure are in good agreement with the results of the two references indicated by the solid 
symbols. The slight variations may be due to differences in finite-difference spacings 
and other approximations in the various numerical procedures. Here the values of the 
rotational spring constant q^ was taken as 10^ to approximate completely constrained 
edge rotation and as 0 for the simply supported edge. 

Meridional Edge Restraint 

Results were also obtained for variations in meridional restraint as well as rota- 
tional restraints. Results were limited to values of ju i 6 since it has been shown that 
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for p > 6, asymmetric buckling governs (ref. 29). To determine the influence of merid- 
ional restraint on the buckling behavior of shallow spherical shells, plots of the buckling 
pressure p cr for various values of the in-plane restraint parameter are shown in 
figure 5 for geometric shell parameters \± of 4, 5, and 6. The calculations were car- 
ried out for a shallow spherical cap with the edge essentially fully restrained against 
rotation (qg = 10^. The plots show that the buckling pressure is essentially insensitive 
to the in-plane restraint at the extreme regions of meridional restraint (q^ < 10 3 and 
qj > 10 7 ). Between these extremes the buckling pressure increases monotonically with 
an increase in in-plane restraint. 

From the above results an estimate can be made of the size of an edge ring of solid 
circular cross section which would be required to approximate a rigid restraint to merid- 
ional displacement. By assuming that q^ = 10 7 and that the ring and shell are made of 
the same material, it is found that a relatively large ring cross-sectional radius of 
approximately 50t is required, where t is the thickness of the shell. 

Rotational Edge Restraint 

Calculations were also made to determine the influence of rotational restraint on 
the buckling behavior of shallow spherical shells. Plots of the buckling pressure p cr 
for various values of the rotational restraint parameter qg are shown in figure 6 for 
geometric parameters p of 4, 5, and 6. The calculations were carried out for a shal- 
low spherical cap with the meridional edge displacement essentially fully restrained 
(q^ = 10 7 ). The plots show that the buckling pressure is essentially insensitive to the 
rotational restraint for the extreme regions of plot (qg < 0.01 and qg > 10 3 ). Between 
these extremes the character of the variation of buckling pressure with a change in the 
rotational restraint parameter is dependent on the value of the geometric shell 
parameter p. 

Figure 6 shows two effects which are of interest. First, that for p = 5 and 6 a 
peak in the p cr curve occurs between the extremes of edge conditions, that is, clamped 
(qg > 10 3 ) and simply supported (qg < 0.01^. Second, that for p = 4 and 5 the buckling 
pressure for the simply supported edge condition is greater than the buckling pressure 
for the clamped edge condition. Although the first effect was unexpected, it can be cor- 
related to some degree with changes in the normal deflection mode shapes. The mode 
shapes just prior to buckling were calculated for several combinations of the parameters. 
The results are shown in figures 7(a), 7(b), and 7(c) for geometric shell parameters p 
of 4, 5, and 6, respectively, and for values of the rotational restraint parameter qg 
associated -with the simply supported and clamped edge conditions. In figures 7(b) 
and 7(c), the modes associated with values of the rotational restraint parameter just 
prior to and just after the peak in p cr are also included. It is seen that the simply 
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supported (qg = o) buckling mode changes between p = 4 and p. = 5 and that the 
clamped » 

With regard to the second effect shown in figure 6, other investigators have noted 
that the spherical shell with a simply supported edge has a higher buckling pressure than 
the same shell with a clamped edge (e.g., ref. 30). The increase in buckling pressure 
was attributed to dynamic effects in reference 30, whereas the results of the present 
investigation show that the increase occurs statically and is associated with the buckling 
mode. Perhaps in this range the buckling pressure for the simply supported edge would 
be lower than for the clamped edge if asymmetric buckling were included in the study. 

The explanation for the peak in the variation in buckling pressure with rotational 
restraint is not apparent. It does appear, however, that the results are not unreasonable 
in view of the variation of buckling pressure with p shown in figure 4. The very sharp 
peak shown in figure 6 for p = 5 is evidently related to the fact that the mode shape for 
the clamped edge has one half-wave whereas the mode shape for the simply supported 
edge has the character of 3 half-waves. 

The size of an edge ring of solid circular cross section required to approximate a 
rigid rotational restraint (q^ > 10^) has been calculated. A ring radius of approximately 
lOt is equivalent to q^ = 103, where the ring and shell are made of the same material. 
This is considerably smaller than the ring radius of approximately 50t required to simu- 
late a rigid restraint to meridional displacement. Thus, in the design of a supporting 
ring to simulate clamped edge conditions, the meridional restraint is the controlling 
factor. 


10 5 ) buckling mode changes between p = 5 and p = 6. 


CONCLUDING REMARKS 

A direct minimization of an object function was used to demonstrate that for non- 
linear static structural problems, stable equilibrium states can be determined numeri- 
cally. The method developed was applied to determine the "top-of-the-knee" buckling 
pressures of shallow spherical shells with various edge restraints. In this case the 
object function was the potential energy of the system. 

Several mathematical programing techniques were evaluated. These techniques 
included steepest-descent, conjugate-gradient, variable -metric and generalized Newton- 
Raphson methods. It was found that the procedure most appropriate for this class of 
problems is a hybrid method composed primarily of a generalized Newton-Raphson pro- 
cedure with selected applications of the conjugate-gradient method. 

This computational procedure was chosen here on the basis that it might be applied 
to some general functional for many classes of problems involving a large number of 
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variables where round-off error might be an important consideration and where maxi- 
mum and inflection points may exist near the desired minimum. 

The procedure has been used to study the influences of a wide range of rotational 
and meridional boundary restraint on the axisymmetric buckling of a shallow spherical 
cap subjected to a uniformly distributed pressure. The variations of buckling pressure 
with changes in in-plane and rotational edge restraint for uniformly loaded shells for 
several values of the geometric shell parameter \i show that for = 5, partial rota- 
tional restraint can yield higher buckling pressures than for either clamped or simply 
supported edges. It was also found that a very large value of meridional in -plane stiff- 
ness was required in order to approach the greater strength corresponding to the com- 
pletely restrained condition. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., December 4, 1969. 
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APPENDIX 


DERIVATIVES OF SHALLOW-SHELL POTENTIAL FUNCTION 


The first and second derivatives of the potential function, equation (8), which are 
required in the Newton-Raphson minimization procedure are found by differentiation of 
the potential with respect to the nodal displacements. 

Gradients 

Differentiation of the linear potential II^ with respect to u^ and w^ leads to 
the following gradient components: 


= e i+(l/2) P i+(l/2)<^i+(l/2) V(l/2)( w l + UW 2) i+(1/2) + 2 W i+d/2) (” i)| 

+ e i-(l/2)Pi-(l/2)^i-(l/2) U i-(l/2)( w l + W 2)._ (1/2) + 2 p 2~|^ W Ml/2)(z)j 
+ e iPi<^i 2w i («J + 4 + 2^ 1 co 2 ) i + 2 ^ 2 + ^i). 2w i(-^) + |f w l(-^r)] 

- 2 fi i} + + ^ + ^ + e i-l p i-l d i-l 2 w Llfe) 


2v 1 


p i-l' 

V 2A 


+ P"^\1a“ + ~^ 2 ~) + 2q 2 w N-l 6 i,N-l + A"( w N 5 i,N ' w N-2 6 i,N-2) 

(i = 1, 2, . . N + 1) (Al) 

= e i+(l/2) p i+(l/2) k i+(l/2) | 2u i + (l/2)(- z} + 2w i+(l/2) (” x)( W l + V ^\ +{l/2 ) 

+ p i+(l/2) Ui+ (l/2) (“ S’) u i+(l/2)(|)| + e i— (l/2) p i— (l/2) k i— (1/2) | 2u l-(l/2)(D 

+ 2w i-(l/2)(l)( w l + W 2)._ (1/2) + p-^ V(1/2)(B + I U i-(l/2)^> + £ i k i p i 2 ^ 


+ 2 pf( w 2 + vw l\ + 


(i = 1, 2, . . N + 1) (A2) 
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where and are in t e g r atinp: factors having values as follows: 


E N+(l/2) 


= 0 


e 4 


i+(l/2) 


€ i = A 


6 1 = 0 

e N+l = 0 

e - A 
e N “ 2 


= A 


(i = 1, 2, . . N) 
(i = 2, 3, . . N - 1) 


and N is the number of stations on the shell. An off-shell station N + 1 is considered 
so that the derivative of (u,w) may be evaluated at the edge. The quantity 5^ is the 


Kronecker delta and is defined as follows: 


6 ij = 


(i*j) 
(i = j) 


In a similar manner contribution of the nonlinear potential to the gradient compo- 
nents is as follows: 


~ k i+(l/2) p i+(l/2) e i+(l/2) Z yMi+(l/2)(" S’) + 2 w 1+(1/2)(~ D W i+(1/2)( W 1 + VU} 2\ 


i+(l/2) 


+ u. 




( W ')i+(l/2) , \ 

+ + w 2). 


+( 1 / 2)1 


+ k i-(l/2) p i-(l/2) e i-(l/2) Z ||( w %-(l/2)(i) + 2w i-(l/2)(i)j 


w. 


i-(l/2)( W l + VW 2\_ 


+ u . 


i-d/2) + " Pii (i/| 


2 V 1 2 A-(1/2)I 


( 1 / 2 ) 


(A3) 
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UST = k i+(l/2) p i+(l/2) e i+(l/2)* <| (w,) i+(l/2) 


V 


A + 2p i+(l/2) 


+ k i-(l/2) p i-(l/2) e i-(l/2) Z f (w,) i-(l/2) 


U 


A 2p i-(l/2) 


(A4) 


where i = 1, 2, . . N + 1. 


Second Derivative of the Potential Function 

The differentiation of the gradient components with respect to nodal displacements 
leads to the following expressions for second derivatives of IIj: 


-f ^Ui/2) ^la/2) U 2 1 2>i i4 1 


(2_ 2v \ 

i+l( A 4 " „ A 3 


V A4 ^i+l A " 


+ d i-l(^A + + 2 ^3 6 i,N-l + T( 5 i,N " 6 i,N-2) 


(A5) 


a2n l -s / 4 2i/ 

d i I “ T7 + 


2d. 


i-(V2) 


9wi awj.i ^ A 4 p . a 3J p 2_ (1/2)A 2 1 " 1 U4 p. iA 3 




2i/ 




(A6) 


2v 


awi aw i+1 


n i - 5-/4 , 2v\ 2d i+(l/2) a / _4_ 

- “ ‘ 1 (a 4 p iA 3 p2 +(1/2)A 2 A 4 Pi+iA 3 y 


(A7) 


8 n l 2d i _ 1 


8Wi 9Wi_ 2 A 4 


(A8) 


9 n l 2d i+1 


^i ^i+2 A 4 


(A9) 


d^u x 

8wi 8uj ^i+(l/2) 


("1 + w 2l „ . ("l + ""2), 2kj , 

A 4/ 1 + k i-(l/2) A ' </) *-^( M 2*-“l) 1 < A10 > 
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a 2 n! 

("i 

- Ir , , 

+ z/o; 0 \ 

2 /i-(l/2) 

(All) 

aui aw^j 

(1/2) 

A 

a 2 ^ 

( W 1 

+ Wn) 

2 Ml/2) 

fA12^ 

9u i ^i+l " 

k i+(l/2) 

A 



a 2 n. 


*1 


2 ~ k i+(l/2) 


2v 


A 2 P i+(1/2) A 


+ k i-(l/2) 


2 + 2v 


_ a2 P i-(1/2) A 


2k, 


+ TT + 2q i 6 i,N-l < A13 > 


920! 2k ; 


814 BUj.j 


i-(l/2) 

A2 


(A14) 


a^IL 2k, 


9Ui 9u i+1 


Ml/2) 

A2 


where 


(A 15) 


k i = k iPiH 
a i = d i p i € i 

Similarly, the second derivatives of the nonlinear potential n n ^ are as follows: 


a 2 !! 


9w 2 


f ~ k i+(l/2) < j 3 ^ w, )i+(l/2)(j) + ^ w i+(l/2)( w l + VOJ 2) i+{1/2) + u i+(l/2) + 


Ml/2) 


p i+(l/2) 


‘ | W i+(l/2)( W l + ^2) i+(l/2) > " k i-(l/2)S 3Z ( w, )i-(l/2)(^2 


A2 


W i (1 /9\f w l + + u! /, /9 x + M 1 / 2 ) 

i“(l/2)l 1 2 4-(l/2) p i-(l/2) 


+ l w !-<i/2)(“i + 1 ’“a) l . (1/2 A 


(A16) 
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S 2 !! 


nl 


9Wi 


= k i-l 
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rX2 (w) Ml/2)-X2 
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nl 
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9wj au£ 
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f 2 t 
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(A20) 
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where i = 1, 2, . . N + 1. 
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